/* ****************************************************************************
 *
 * Copyright (c) Microsoft Corporation. 
 *
 * This source code is subject to terms and conditions of the Microsoft Public License. A 
 * copy of the license can be found in the License.html file at the root of this distribution. If 
 * you cannot locate the  Microsoft Public License, please send an email to 
 * dlr@microsoft.com. By using this source code in any fashion, you are agreeing to be bound 
 * by the terms of the Microsoft Public License.
 *
 * You must not remove this notice, or any other, from this software.
 *
 * THIS CODE AND INFORMATION ARE PROVIDED "AS IS" WITHOUT WARRANTY OF ANY
 * KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A
 * PARTICULAR PURPOSE.
 * ***************************************************************************/


#if !defined(AFX_CONTOUR_H__45A84754_43CD_414A_80A0_9211F954B0DA__INCLUDED_)
#define AFX_CONTOUR_H__45A84754_43CD_414A_80A0_9211F954B0DA__INCLUDED_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000

#include "Field.h"


//! Contouring algorithm 
namespace contour
{

/*!  \brief A contouring class

  TContour implements Contour plot algorithm descrided in 
		IMPLEMENTATION OF
		AN IMPROVED CONTOUR
		PLOTTING ALGORITHM
		BY 

		MICHAEL JOSEPH ARAMINI 

		B.S., Stevens Institute of Technology, 1980 
 See http://www.ultranet.com/~aramini/thesis.html

 Ported to C++ by Jonathan de Halleux.

 Using TContour :

 TContour is not directly usable. The user has to 
	-1. derive the function ExportLine that is 
		supposed to draw/store the segment of the contour
	-2. Attach a field object to it

\ingroup ContourLibGroup

*/
template<class T>
class TContour  
{
public:
	typedef std::vector<T> PlaneContainer;
	struct SFnStr 
	{
	    T m_dFnVal;
	    short m_sLeftLen;
		short m_sRightLen;
		short m_sTopLen;
		short m_sBotLen;
	};

	//! \name Constructors
	//@{
	//! Default constructor
	TContour();
	virtual ~TContour();
	//@}

	//! \name Main Methods
	//@{
	//! Initialize memory. Called in Generate
	virtual void InitMemory();
	//! Clean work arrays
	virtual void CleanMemory();

	/*! \brief Generates contour

	Before calling this functions you must
		-1. derive the function ExportLine that is 
			supposed to draw/store the segment of the contour
		-2. Set the function draw contour of. (using  SetFieldFn
			The function must be declared as follows
			double (*myF)(double x , double y); 
	*/
	virtual void Generate();
	//@}

	//! \name Grids
	//@{
	//! Sets the pointer to the F(x,y) funtion
	void SetField(TField<T>* _pField);
	//! Returns pointer to the current field function
	TField<T>* GetField()					{	return m_pField;};
	//@}

	//! \name Grid Methods
	//@{
	//! Set the dimension of the primary grid
	void SetFirstGrid(int iRow, int iCol);
	//! Set the dimension of the base grid
	void SetSecondaryGrid(int iRow, int iCol);
	/*! \brief Sets the number of isocurve and their values
	 pPlanes must be allocated on the heap. It will be deleted by TContour.
	 */
	//! Retrieve number of columns of primary grid
	int GetColFir() const		{	return m_iColFir;};
	//! Retrieve number of rows of primary grid
	int GetRowFir() const		{	return m_iRowFir;};
	//! Retrieve number of columns of secondary grid
	int GetColSec() const		{	return m_iColSec;};
	//! Retrieve number of rows of secondary grid
	int GetRowSec() const		{	return m_iRowSec;};
	//@}

	//! \name Plane methods
	//@{
	void SetPlanes(const PlaneContainer vPlanes){	CleanMemory();m_vPlanes=vPlanes;};
	//! Return current number of planes
	size_t GetPlanesSize() const		{ return m_vPlanes.size();};
	//! Returns planes
	const PlaneContainer& GetPlanes()	const {	return m_vPlanes;};
	//! Returns i-th plane
	double GetPlane(UINT i)	const	{	ASSERT(i>=0); ASSERT(i<m_vPlanes.size()); return m_vPlanes[i];};
	//@}

	//! \name Approximated point coordinates
	//@{
	//! For an indexed point i on the sec. grid, returns x(i)
	T GetXi(int i) const	
	{	
		const TField<T>::SLimits& limits = m_pField->GetLimits(); 
		return (T)(limits.left+(i%(m_iColSec+1))*(limits.GetWidth())/(T)( m_iColSec ));
	};
	//! For an indexed point i on the fir. grid, returns y(i)
	T GetYi(int i) const	
	{
		const TField<T>::SLimits& limits = m_pField->GetLimits(); 		
		UINT index=i/(m_iColSec+1);
		return (T)(limits.bottom+index*limits.GetHeight()/(T)(m_iRowSec));
	};
	//@}

protected:
	//! \name Attributes
	//@{
	//! value of contour planes
	PlaneContainer m_vPlanes;			
	//! primary	grid, number of columns
	int m_iColFir;	
	//! primary	grid, number of rows
	int m_iRowFir;
	//! secondary grid, number of columns
	int m_iColSec;		
	//! secondary grid, number of rows
	int m_iRowSec;				
	T m_dDx;
	T m_dDy;
	//! pointer to F(x,y) function
	TField<T>* m_pField;
	//@}

	//! Protected function
	virtual void ExportLine(UINT iPlane,int x1, int y1, int x2, int y2) = 0; // plots a line from (x1,y1) to (x2,y2)

	// Work functions and variables
	//! pointer to mesh parts
	SFnStr** m_ppFnData;
	SFnStr* FnctData(int i,int j)			{	return (m_ppFnData[i]+j);};
	T Field(int x, int y);	 /* evaluate funct if we must,	*/
	void Cntr1(int x1, int x2, int y1, int y2);
	void Pass2(int x1, int x2, int y1, int y2);	  /* draws the contour lines */
};

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

template <class T>
TContour<T>::TContour()
: m_vPlanes(1,0)
{
	m_iColFir=m_iRowFir=32;
	m_iColSec=m_iRowSec=256;
	m_dDx=m_dDy=0;
	m_pField=NULL;
	m_ppFnData=NULL;
}

template <class T>
TContour<T>::~TContour()
{
	CleanMemory();
}

template <class T>
void TContour<T>::InitMemory()
{
	if (!m_ppFnData)
	{
		m_ppFnData=new SFnStr*[m_iColSec+1];
		for (int i=0;i<m_iColSec+1;i++)
		{
			m_ppFnData[i]=NULL;
		}
	}
}

template <class T>
void TContour<T>::CleanMemory()
{
	if (m_ppFnData)
	{
		int i;
		for (i=0;i<m_iColSec+1;i++)
		{
			if (m_ppFnData[i])
				delete[] (m_ppFnData[i]);
		}
		delete[] m_ppFnData;
		m_ppFnData=NULL;
	}
}

template <class T>
void TContour<T>::Generate()
{
	if (!m_pField)
		return;

	int i, j;
	int x3, x4, y3, y4, x, y, oldx3, xlow;
	const int cols=m_iColSec+1;
	const int rows=m_iRowSec+1;
	T xoff,yoff;

	// Initialize memroy if needed
	InitMemory();

	const TField<T>::SLimits& limits = m_pField->GetLimits();	

	m_dDx = limits.GetWidth()/(T)(m_iColSec);
	xoff = limits.left;
	m_dDy = limits.GetHeight()/(T)(m_iRowSec);
	yoff = limits.bottom;

	xlow = 0;
	oldx3 = 0;
	x3 = (cols-1)/m_iRowFir;
	x4 = ( 2*(cols-1) )/m_iRowFir;
	for (x = oldx3; x <= x4; x++) 
	{	  /* allocate new columns needed
		*/
		if (x >= cols)
			break;
		if (m_ppFnData[x]==NULL)
			m_ppFnData[x] = new SFnStr[rows];

		for (y = 0; y < rows; y++)
			FnctData(x,y)->m_sTopLen = -1;
	}

	y4 = 0;
	for (j = 0; j < m_iColFir; j++) 
	{
		y3 = y4;
		y4 = ((j+1)*(rows-1))/m_iColFir;
		Cntr1(oldx3, x3, y3, y4);
	}

	for (i = 1; i < m_iRowFir; i++) 
	{
		y4 = 0;
		for (j = 0; j < m_iColFir; j++) 
		{
			y3 = y4;
			y4 = ((j+1)*(rows-1))/m_iColFir;
			Cntr1(x3, x4, y3, y4);
		}

		y4 = 0;
		for (j = 0; j < m_iColFir; j++) 
		{
			y3 = y4;
			y4 = ((j+1)*(rows-1))/m_iColFir;
			Pass2(oldx3,x3,y3,y4);
		}

		if (i < (m_iRowFir-1)) 
		{	 /* re-use columns no longer needed */
			oldx3 = x3;
			x3 = x4;
			x4 = ((i+2)*(cols-1))/m_iRowFir;
			for (x = x3+1; x <= x4; x++) 
			{
				if (xlow < oldx3) 
				{
					if (m_ppFnData[x])
						delete[] m_ppFnData[x];
					m_ppFnData[x] = m_ppFnData[xlow];
					m_ppFnData[ xlow++ ] = NULL;
				} 
				else
					if (m_ppFnData[x]==NULL)
						m_ppFnData[x] = new SFnStr[rows];

				for (y = 0; y < rows; y++)
					FnctData(x,y)->m_sTopLen = -1;
			}
		}
	}

	y4 = 0;
	for (j = 0; j < m_iColFir; j++) 
	{
		y3 = y4;
		y4 = ((j+1)*(rows-1))/m_iColFir;
		Pass2(x3,x4,y3,y4);
	}
}

template <class T>
void TContour<T>::Cntr1(int x1, int x2, int y1, int y2)
{
	T f11, f12, f21, f22, f33;
	int x3, y3, i, j;
	
	if ((x1 == x2) || (y1 == y2))	/* if not a real cell, punt */
		return;
	f11 = Field(x1, y1);
	f12 = Field(x1, y2);
	f21 = Field(x2, y1);
	f22 = Field(x2, y2);
	if ((x2 > x1+1) || (y2 > y1+1)) {	/* is cell divisible? */
		x3 = (x1+x2)/2;
		y3 = (y1+y2)/2;
		f33 = Field(x3, y3);
		i = j = 0;
		if (f33 < f11) i++; else if (f33 > f11) j++;
		if (f33 < f12) i++; else if (f33 > f12) j++;
		if (f33 < f21) i++; else if (f33 > f21) j++;
		if (f33 < f22) i++; else if (f33 > f22) j++;
		if ((i > 2) || (j > 2)) /* should we divide cell? */
		{	
			/* subdivide cell */
			Cntr1(x1, x3, y1, y3);
			Cntr1(x3, x2, y1, y3);
			Cntr1(x1, x3, y3, y2);
			Cntr1(x3, x2, y3, y2);
			return;
		}
	}
	/* install cell in array */
	FnctData(x1,y2)->m_sBotLen = FnctData(x1,y1)->m_sTopLen = x2-x1;
	FnctData(x2,y1)->m_sLeftLen = FnctData(x1,y1)->m_sRightLen = y2-y1;
}

template <class T>
void TContour<T>::Pass2(int x1, int x2, int y1, int y2)
{
	int left, right, top, bot,old, iNew, x3, y3;
	UINT i,j;
	T yy0, yy1, xx0, xx1, xx3, yy3;
	T v, f11, f12, f21, f22, f33, fold, fnew, f;
	_ASSERT(m_pField);
	const TField<T>::SLimits& limits = m_pField->GetLimits();	

	T xoff=limits.left;
	T yoff=limits.bottom;
	
	if ((x1 == x2) || (y1 == y2))	/* if not a real cell, punt */
		return;
	f11 = FnctData(x1,y1)->m_dFnVal;
	f12 = FnctData(x1,y2)->m_dFnVal;
	f21 = FnctData(x2,y1)->m_dFnVal;
	f22 = FnctData(x2,y2)->m_dFnVal;
	if ((x2 > x1+1) || (y2 > y1+1)) /* is cell divisible? */
	{	
		x3 = (x1+x2)/2;
		y3 = (y1+y2)/2;
		f33 = FnctData(x3, y3)->m_dFnVal;
		i = j = 0;
		if (f33 < f11) i++; else if (f33 > f11) j++;
		if (f33 < f12) i++; else if (f33 > f12) j++;
		if (f33 < f21) i++; else if (f33 > f21) j++;
		if (f33 < f22) i++; else if (f33 > f22) j++;
		if ((i > 2) || (j > 2)) /* should we divide cell? */ 
		{	
			/* subdivide cell */
			Pass2(x1, x3, y1, y3);
			Pass2(x3, x2, y1, y3);
			Pass2(x1, x3, y3, y2);
			Pass2(x3, x2, y3, y2);
			return;
		}
	}

	for (i = 0; i < m_vPlanes.size(); i++) 
	{
		v = m_vPlanes[i];
		j = 0;
		if (f21 > v) j++;
		if (f11 > v) j |= 2;
		if (f22 > v) j |= 4;
		if (f12 > v) j |= 010;
		if ((f11 > v) ^ (f12 > v)) 
		{
			if ((FnctData(x1,y1)->m_sLeftLen != 0) &&
				(FnctData(x1,y1)->m_sLeftLen < FnctData(x1,y1)->m_sRightLen)) 
			{
				old = y1;
				fold = f11;
				while (1) 
				{
					iNew = old+FnctData(x1,old)->m_sLeftLen;
					fnew = FnctData(x1,iNew)->m_dFnVal;
					if ((fnew > v) ^ (fold > v))
						break;
					old = iNew;
					fold = fnew;
				}
				yy0 = ((old-y1)+(iNew-old)*(v-fold)/(fnew-fold))/(y2-y1);
			} 
			else
				yy0 = (v-f11)/(f12-f11);

			left = (int)(y1+(y2-y1)*yy0+0.5);
		}
		if ((f21 > v) ^ (f22 > v)) 
		{
			if ((FnctData(x2,y1)->m_sRightLen != 0) &&
				(FnctData(x2,y1)->m_sRightLen < FnctData(x2,y1)->m_sLeftLen)) 
			{
				old = y1;
				fold = f21;
				while (1) 
				{
					iNew = old+FnctData(x2,old)->m_sRightLen;
					fnew = FnctData(x2,iNew)->m_dFnVal;
					if ((fnew > v) ^ (fold > v))
						break;
					old = iNew;
					fold = fnew;
				}
				yy1 = ((old-y1)+(iNew-old)*(v-fold)/(fnew-fold))/(y2-y1);
			} 
			else
				yy1 = (v-f21)/(f22-f21);

			right = (int)(y1+(y2-y1)*yy1+0.5);
		}
		if ((f21 > v) ^ (f11 > v)) 
		{
			if ((FnctData(x1,y1)->m_sBotLen != 0) &&
				(FnctData(x1,y1)->m_sBotLen < FnctData(x1,y1)->m_sTopLen)) {
				old = x1;
				fold = f11;
				while (1) {
					iNew = old+FnctData(old,y1)->m_sBotLen;
					fnew = FnctData(iNew,y1)->m_dFnVal;
					if ((fnew > v) ^ (fold > v))
						break;
					old = iNew;
					fold = fnew;
				}
				xx0 = ((old-x1)+(iNew-old)*(v-fold)/(fnew-fold))/(x2-x1);
			} 
			else
				xx0 = (v-f11)/(f21-f11);

			bot = (int)(x1+(x2-x1)*xx0+0.5);
		}
		if ((f22 > v) ^ (f12 > v)) 
		{
			if ((FnctData(x1,y2)->m_sTopLen != 0) &&
				(FnctData(x1,y2)->m_sTopLen < FnctData(x1,y2)->m_sBotLen)) {
				old = x1;
				fold = f12;
				while (1) {
					iNew = old+FnctData(old,y2)->m_sTopLen;
					fnew = FnctData(iNew,y2)->m_dFnVal;
					if ((fnew > v) ^ (fold > v))
						break;
					old = iNew;
					fold = fnew;
				}
				xx1 = ((old-x1)+(iNew-old)*(v-fold)/(fnew-fold))/(x2-x1);
			} 
			else
				xx1 = (v-f12)/(f22-f12);

			top = (int)(x1+(x2-x1)*xx1+0.5);
		}

		switch (j) 
		{
			case 7:
			case 010:
				ExportLine(i,x1,left,top,y2);
				break;
			case 5:
			case 012:
				ExportLine(i,bot,y1,top,y2);
				break;
			case 2:
			case 015:
				ExportLine(i,x1,left,bot,y1);
			break;
		case 4:
		case 013:
			ExportLine(i,top,y2,x2,right);
			break;
		case 3:
		case 014:
			ExportLine(i,x1,left,x2,right);
			break;
		case 1:
		case 016:
			ExportLine(i,bot,y1,x2,right);
			break;
		case 0:
		case 017:
			break;
		case 6:
		case 011:
			yy3 = (T)((xx0*(yy1-yy0)+yy0)/(1.0-(xx1-xx0)*(yy1-yy0)));
			xx3 = yy3*(xx1-xx0)+xx0;
			xx3 = x1+xx3*(x2-x1);
			yy3 = y1+yy3*(y2-y1);
			xx3 = xoff+xx3*m_dDx;
			yy3 = yoff+yy3*m_dDy;

			_ASSERT(m_pField);
			f = m_pField->Eval(xx3, yy3);
			if (f == v) {
				ExportLine(i,bot,y1,top,y2);
				ExportLine(i,x1,left,x2,right);
			} else
				if (((f > v) && (f22 > v)) || ((f < v) && (f22 < v))) {
					ExportLine(i,x1,left,top,y2);
					ExportLine(i,bot,y1,x2,right);
				} else {
					ExportLine(i,x1,left,bot,y1);
					ExportLine(i,top,y2,x2,right);
				}
		}
	}
}

template <class T>
T TContour<T>::Field(int x, int y)	 /* evaluate funct if we must,	*/
{
	T x1, y1;
	_ASSERT(m_pField);
	const TField<T>::SLimits& limits = m_pField->GetLimits();	
	
	if (FnctData(x,y)->m_sTopLen != -1)  /* is it already in the array */
		return(FnctData(x,y)->m_dFnVal);

	/* not in the array, create new array element */
	x1 = limits.left+m_dDx*x;
	y1 = limits.bottom+m_dDy*y;
	FnctData(x,y)->m_sTopLen = 0;
	FnctData(x,y)->m_sBotLen = 0;
	FnctData(x,y)->m_sRightLen = 0;
	FnctData(x,y)->m_sLeftLen = 0;

	return (FnctData(x,y)->m_dFnVal = m_pField->Eval(x1, y1));
}

template <class T>
void TContour<T>::SetFirstGrid(int iRow, int iCol)
{
	m_iColFir=__max(iCol,2);
	m_iRowFir=__max(iRow,2);
}

template <class T>
void TContour<T>::SetSecondaryGrid(int iRow, int iCol)
{
	// cleaning work matrices if allocated
	CleanMemory();

	m_iColSec=__max(iCol,2);
	m_iRowSec=__max(iRow,2);
}

template <class T>
void TContour<T>::SetField(TField<T>* _pField)
{	
	m_pField=_pField; 
};

};

#endif // !defined(AFX_CONTOUR_H__45A84754_43CD_414A_80A0_9211F954B0DA__INCLUDED_)
